/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file

    Implementation of the Gadgetheader routines. */

//$Id$

#include "gadgetquant.h"
#include <iostream>
#include "biniostream.h"

using namespace std;

/** Loads the header from the specified Gadget type 1 snapshot
    file. Since some header fields have different meaning depending on
    whether it's a "Springel" or "TJ" type file, a heuristic is used
    to try to figure this out, but it is not robust. */
void mcrx::Gadgetheader::loadfrheader(const char snapfn[255])
{
  int k;
  binifstream gsnap;
  gsnap.open(snapfn);
  gsnap>>k;
  bool byte_swap=false;
  if (k == 65536) {
    gsnap.swap(); 
    byte_swap=true;
  }

  int ntot=0;
  for(k=0; k<=5; k++)   {
      gsnap >> npart[k];
      ntot+=npart[k];
    }

  for(k=0; k<=5; k++)
    {
      gsnap >> mass[k];
    }
  gsnap >> time;
  gsnap >> redshift;
  gsnap >> flag_sfr;
  gsnap >> flag_feedback;

  for(k=0; k<=5; k++)
    {
      gsnap >> npartTotal[k];
    }

  gsnap >> flag_cooling;
  gsnap >> num_files;
  gsnap >> BoxSize;
  gsnap >> Omega0;
  gsnap >> OmegaLambda;
  gsnap >> HubbleParam;

  // conditionally read next to fields
  int temp1,temp2;
  gsnap >> temp1 >> temp2;

  // check for Springel file.  the heuristic used is to check the
  // position of the metallicity and potential fields
  {
    // we now proceed on the ASSUMPTION that it is a springel file
    // and see if it makes sense
    int age=temp1,metals=temp2;
    int nmass=0;
    for(int k=0; k<=5; k++)
      nmass+= (mass[k]==0)?npart[k]:0;
    binifstream seek(snapfn);
    if(byte_swap)
      seek.swap();
    const int offset=264+7*4*ntot+nmass*4+
      4*npart[0]*(3+2*flag_cooling+flag_sfr)+
      4*npart[4]*age+
      8*(7+2*flag_cooling+flag_sfr+int(age&&(npart[4]>0)) );

    seek.skip(offset);
    // if this is the springel file, we should not be looking at
    // metallicity or potential, depending on whether metallicity is
    // set
    int test;
    bool gotit=false;
    if(metals) {
      seek >> test;
      int next=4*(npart[0]+npart[4]);
      if(test == next) {
	//  Field looks like springel, continue
	seek.skip(next+4);
      }
      else {
	// metallicity does not look like Springel
	is_springel_file=false;
	gotit=true;
      }
    }
    // now we should have pot
    if (!gotit) {
      seek >> test;
      int next=4*ntot;
      if(test == next) {
	//  Field looks like springel, continue
	seek.skip(next+4);
      }
      else {
	// does not look like Springel
	is_springel_file=false;
	gotit=true;
      }
    }
    // now there are mysterious short fields for the springel file...
    test=0;
    if (!gotit) {
      seek >> test;
      if((test > 100) && seek.good()) {
	// does not look like Springel
	is_springel_file=false;
	gotit=true;
      }
      else
	// consistent with Springel
	is_springel_file=true;
    }
  }

  if(is_springel_file)
    cerr << "GADGET file has Springel format" << endl;
  else
    cerr << "GADGET file has TJ format" << endl;

  if(is_springel_file) {
    flag_multiphase=0;
    flag_stellarage=temp1;
    flag_sfrhistogram=0;
    flag_stargens=1;
    flag_snapaspot=1;
    flag_metals=temp2;
    flag_energydetail=0;
    flag_parentID=0;
    flag_starorig=0;
  }
  else {
    flag_multiphase=temp1;
    flag_stellarage=temp2;
    gsnap >> flag_sfrhistogram;
    gsnap >> flag_stargens;
    gsnap >> flag_snapaspot;
    gsnap >> flag_metals;
    gsnap >> flag_energydetail;
    gsnap >> flag_parentID >> flag_starorig;
  }
}

void mcrx::Gadgetheader::printheader() const
{
  int k;
  cout << "SNAPSHOT HEADER QUANTITIES\n";
  cout << "Number of particles:\n";
  for(k=0; k<=5; k++)   
    {
      cout <<  npart[k] << endl;
    }
  cout << "Particle masses\n";
  for(k=0; k<=5; k++)
    {
      cout << mass[k] << endl;
    }
  cout << "Time: " << time << endl;
  cout << "z: " << redshift << endl;
  cout << "flag_sfr: " << flag_sfr << endl;
  cout << "flag_feedback: " << flag_feedback << endl;
  cout << "Total number of particles\n";
  for(k=0; k<=5; k++)
    {
      cout << npartTotal[k] << endl;
    }
  cout << "flag_cooling: " << flag_cooling << endl;
  cout << "num_files: " << num_files << endl;
  cout << "BoxSize: " << BoxSize << endl;
  cout << "Omega0: " << Omega0 << endl;
  cout << "OmegaLambda: " << OmegaLambda << endl;
  cout << "HubbleParam: " << HubbleParam << endl;
  cout << "flag_multiphase: " << flag_multiphase << endl;
  cout << "flag_stellarage: " << flag_stellarage << endl;
  cout << "flag_sfrhistogram: " << flag_sfrhistogram << endl; 
  cout << "flag_stargens: " << flag_stargens << endl; 
  cout << "flag_snapaspot: " << flag_snapaspot << endl; 
  cout << "flag_metals: " << flag_metals << endl; 
  cout << "flag_energydetail: " <<flag_energydetail << endl; 
}










